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ABSTRACT 


We present the results of a search for correlations between X-ray 
flux and angular acceleration for the accreting binary pulsar Vela X-l, 
based on data obtained with the Hakucho satellite during the interval 
1982-1984. In undertaking this correlation analysis, we found it necessary 
to modify the usual statistical method to deal with conditions imposed 
by generally unavoidable satellite observing constraints, most notably 
a mismatch in sampling between the two variables. Our results are 
suggestive of a correlation between flux and the absolute value of the 
angular acceleration, at a significance level of 96 percent. We discuss 
the implications of the methods and results presented here for future 
observations and analysis. 

Subject headings: stars: neutron stars: rotation X-rays: binaries 

pulsars numerical methods 



( 3 ) 

I. INTRODUCTION 


Variations in the rotation rate of X-ray pulsars are an inevitable 
consequence of matter accretion. Theoretical studies indicate that a 
description of accretion torques must take into account not only the 
mass accretion rate, but the accretion flow pattern as well as the 
role of the neutron star magnetosphere (Lamb 1989, and references 
therein) . Observations of the variation in rotation rate for several 
accretion-powered pulsars have provided significant constraints on these 
models (Boynton et al. 1984; Nagase et al. 1984; Deeter et al. 1989). 
However, a more detailed probing of accretion physics is possible by 
simultaneously observing both the neutron star angular acceleration, a, 
and the X-ray flux, F, the two principal manifestations of the accretion 
process. According to preliminary theoretical studies and simulations 
reported by Lamb (1985), the character of the temporal correlations between 
F and a is distinctly dependent on whether the accretion flow is radial or 
disk-mediated, and whether the system can be described as a fast or a slow 
rotator . 

Despite the attractiveness of investigating correlations between F 
and a as a tool for understanding accretion flows, only one observational 
study has previously been reported, for the transient X-ray source 
EXO 2030+375 (Parmar et al . 1988) . This is partly because of difficulty 
in obtaining adequate observations. For most pulsars, the angular 

acceleration is much smaller than its uncertainty for intervals shorter 
than a few days and it therefore takes a demanding observing program to 
obtain even a small number of independent (F, a) pairs. But even with an 
adequate data base, correctly accounting for the fundamental disparity in 
the way that F and a are sampled requires the development of an extention 
of the customary statistical methods for analyzing and interpreting the 
correlations between these variables. 

This paper consequently has two purposes: to explicate and modify 

the customary mathematical and statistical tools necessary to undertake a 
determination of the correlation between flux and angular acceleration, 
and then to apply the resulting procedures to the Hakucho observations 


of the accretion-powered pulsar Vela X-l . The result of this analysis 
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suggests a possible correlation between F and a , but this result is by 
no means definitive. However, as a consequence of developing a method of 
analysis, we are able to indicate conditions for future observations of 
other sources that should provide a more favorable opportunity to study 
accretion physics through F-c* correlations. 

In § II we briefly review a simple statistical model based on the 
physics relating X-ray flux and neutron star angular acceleration . In 
§ III we outline the method of analysis used in the correlation study, 
leaving details for a series of Appendices. In § IV we discuss available 
data sets and our choice for the analysis in this paper. In § V we present 
the results of our study of the Hakucho data on Vela X-l. In § VI we 
discuss the implications of these results for future observations and 
analysis . 


II. FOUNDATIONS FOR A STATISTICAL MODEL 

In the conventional interpretation, an X-ray pulsar is a rotating 
neutron star whose luminosity arises from mass accretion. Because the 
accreted material carries angular momentum, there will be an accretion 
torque and consequently changes in the rotation rate of the neutron star. 
There is sufficient theoretical understanding of the resulting correlation 
between mass accretion rate and the accretion torque to distinguish between 
various types of accretion flows (Lamb 1985; Lamb et al. 1989). 

What makes F and a the appropriate choice of observables for studying 
the physics of accretion? The physical models proposed by Ghosh and 
Lamb (1978, 1979a, b) can all be expressed as relations between the mass 
accretion rate M and the torque N acting on the neutron star. Under 
rather mild assumptions, these can be restated as relations between the 
more directly observable variables F and a. In particular, the star 
must be sufficiently rigid so that internal torques are unimportant and 
thus a = N/I s , where I s is the moment of inertia of the neutron star. 

On the other hand, the X-ray flux must be a good indicator of the mass 
accretion rate (no beam wandering, no variable absorption, no delayed 
nuclear burning of accreted material). Under these conditions, the a-F 
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relation is identical up to a multiplicative constant to the N-M relation. 
This relationship can be established for a particular source provided the 
accretion rate is sufficiently variable to allow sampling of a reasonable 
range of both F and a. 

In fact, we observe the angular acceleration of the neutron star 
rotation indirectly, identifying it with the time derivative of the X-ray 
pulse frequency. This identification may be comprimised by pulse shape 
noise and by wandering of the X-ray beam pattern. We discuss both of 
these issues later in this section, justifying the naive point of view that 
pulse frequency is identical to rotation frequency, and that the measured 
frequency differs from both only because of uncertainties introduced by 
pulse shape noise. This convention allows us to keep the notation simple, 
using a(t) to denote the time history of the angular acceleration of the 
neutron star rotation, and a (or a similar modification of a) for the 
measured derivative of the pulse frequency averaged over some specified 
interval. Thus, in the remainder of this paper, a will formally represent 
the derivative of the pulse frequency (the observed variable) , which we 
interpret whenever necessary as the angular acceleration of the neutron 
star rotation . 

Except for the transient source EXO 2030+375 (Parmar et al. 1988) , 
direct observational evidence for a relationship between F and a is 
presently rather weak for any X-ray pulsar, but there is substantial 
indirect evidence that these correlations could be present in Vela 
X-l . For instance, both flux and angular acceleration exhibit large 
fluctuations, and the noise properties of the two variables are strikingly 
similar. This similarity is virtually a sine qua non for the 
existence of correlations, since any two correlated variables should 
have nearly identical noise properties. In the following discussion of 
the fluctuations in F and a, we concentrate on comparing the similaries in 
their power density spectra. 

Vela X-l exhibits large variations in angular acceleration with 
spindown and spinup observed on a wide range of timescales (Nagase 
et al. 1984; Boynton et al. 1984; Deeter et al. 1989) . The power density 
spectrum of a has been previously discussed by Boynton et al. (1984) and 
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Deeter et al. (1987, 1989), who show that it is nearly flat for analysis 

frequencies between 0.0003 and 0.2 d — ^ , and follows approximately an f^ 
power law for higher frequencies. The interpretion of the spectrum 

of a as a superposition of two distinct noise processes having spectra 
described by two different power laws (Deeter et al. 1989) is a crucial 
ingredient in a comparision with the spectrum for F . At analysis 
frequencies below 0.2 d — ^ the power density spectrum for a is dominated by 
variations in the rotation rate of the neutron star, and this portion of 
the spectrum is flat (f® power law) within the uncertainty in determining 
the power-law exponent. For analysis frequencies above the crossover 
frequency 0.2 d~ a separate process with an f^ power law dominates the 
spectrum. We understand this second component as arising from fluctuations 
in pulse shape . inducing fluctuations in the measured pulse phase, 
frequency, and angular acceleration. The division of the power density 
spectrum of the variations in a into two regimes is particularly appealing 
because of the simplicity of explaining the two different power laws in the 
low and high frequency regimes in terms of distinct processes. 

There are, however, possible ambiguities in this simple interpretation 
of the spectrum, particularly at analysis frequencies in the vicinity of 
the transition between the two power laws. For instance, van der Klis 
and Bonnet-Bidaud (1984) have suggested that part of the variations at 
lower analysis frequencies might be due to changes in the direction of the 
pulsar beam pattern rather than true variations in the rotation rate of the 
neutron star. This explanation can be ruled out at analysis frequencies 
lower than 0.01 d — ^ , for which the inferred phase excursions become larger 
than 7r radians (Deeter et al. 1989) . It is therefore almost certain that 
the fluctuations at the lowest frequencies are due to variations in the 
rotation rate of the crust of the neutron star. At frequencies between 
0.01 and 0.2 d — ^ , however, a substantial fraction of the power density 
could be due to beam wandering. In this region the spectrum of a is still 
nearly flat, and any power density arising from beam wandering must be 
somewhat compensated by a corresponding decrease in the power density 
from fluctuations in the rotation rate. Even though there is a degree of 
articif iality in proposing a third component in the fluctuations of the 
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measured angular acceleration, further evidence is required to rule out 
this case altogether . 

To strengthen the argument against a third component in the spectrum 
of the fluctuations in a, we consider what observable effect beam 
wandering would have on the pulse shape . In the simplest case of a beam 
with a fairly narrow, two-dimensional gaussian profile, no change in 
pulse shape would result from wandering in azimuth and only a change in 
amplitude would result from wandering in elevation, where azimuth and 
elevation refer to a fixed coordinate system centered on the neutron 
star. This simple behavior clearly does not apply for pulsations with 
complex waveforms, such as seen in Vela X-l, for which it is difficult 
to have beam wandering without a concomitant change in pulse shape. In 
particular, changes in pulse shape entail changes in the relative phases of 
the pulse harmonics, and we can directly compare the fluctuations in these 
harmonic phase differences to the fluctuations in overall pulse phase. 

Thus, in Figure 1, we show the variation with time of the pulse phase of 
Vela X-l for a subset of the Hakucho data, together with the phases of 
the individual harmonic components relative to a smooth trend obtained by 
fitting a quadratic function in time to the variation in pulse phase. In 
this figure, the first harmonic has been omitted, because its amplitude is 
close to zero, and the phase of the third harmonic is relatively noisy as 
its amplitude is small. The remaining harmonic phases track each other 
quite well for time scales longer than about 2 days (corresponding to 
analysis frequencies lower than 0.08 d~*), indicating no appreciable 
change in pulse shape despite large changes in pulse phase over this same 
time interval. ^ This largely closes the frequency gap within which beam 


^ There are hints of very small secular variations in the differential 
phases, but these are at least an order of magnitude smaller than the 
change in overall pulse phase. 

wandering can contribute to the power density spectrum of the phase (and 
angular acceleration) fluctuations. 

We turn now to variations in the X-ray flux, for which fluctuations of 
a factor of five are obseved on a timescale of hours (Nagase et al. 1984). 
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Using the method presented in Appendix A, we constructed a power density 
spectrum of the flux for analysis frequencies between 0.001 and 4 d~ 
based on data obtained with Hakucho . This spectrum, shown in Figure 2, 
is virtually flat for analysis frequencies between 0.001 and 0.3 d“^. 

In contrast to a, all variations in the flux, such as those due to the 
finite number of photons detected, are likely to have a flat power density 
spectrum. Thus, the characteristic crossover frequency observed in the 
spectrum of the fluctuations in a is not present in the corresponding 
spectrum for F. The decrease in power density at the higher frequencies 
is in agreement with the nature of the variations seen in the time history, 
wherein the fluctuations are relatively smooth on time scales smaller than 
one or two hours (Nagase et al . 1984) . 

The power density spectra of the fluctuations in F and a are very 
similar within their common range of analysis frequencies from 0.001 to 
0.2 d“^. Above the crossover frequency 0.2 d~“^, where fluctuations in 
the rotation rate of the neutron star are obscured by pulse shape noise, 
we cannot directly compare the two spectra. Thus we cannot discern 
whether the spectrum of the rotational angular acceleration decreases above 
0.3 d”^ as does the spectrum of F. 

The similarity of the two power density spectra is a precondition for 
a strong correlation between F and a, but of course does not by itself 
establish the existence of a correlation. To determine whether or not 
a correlation is present in a particular data set requires the use of a 
statistical correlation analysis technique. Even so, the power density 
spectra provide information that is helpful in undertaking and interpreting 
the correlation analysis. For instance, the drop in flux power density 
at frequencies above 0.3 d“* , indicating that the flux variations are 
smooth on timescales shorter than a few hours, but that they undergo large 
excursion on longer timescales. On the other hand, determinations of a 
on single intervals shorter than 3 or 4 days are dominated by noise from 
variations in pulse shape, and are therefore useless for a correlation 
diagram. To get correlation pairs (F, a) with a meaningful estimate for a, 
we for shorter intervals we must average over many separate excursions in 
the mass accretion rate. 
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Averaging does not wash out any approximately linear relationship 
between the two variables, but the magnitude of the variations will 
likely decrease as the averaging time scale is lengthened. However, 
averaging over many independent “events” ( separate excursions in the 
mass accretion rate) suggests a statistical approach to analyzing the flux 
and angular acceleration for correlations. In § III we develop suitable 
technique, showing how to modify a common statistical treatment to deal 
with the several complications inherent to these data. 

III. METHODOLOGY 

a) Statistic for the Correlation Coefficient 

For simple accretion-torque models, variations in X-ray flux and 
in the rotational acceleration are both correlated with variations in 
mass accretion, and hence with each other. We observe that variations 
in F and a separately have the character of white noise, with a possible 
high-frequency cut-off set by the duration of an individual accretion 
event. Averages of F and a over an interval whose length T is larger than 
the duration of accretion events have variances which scale inversely with 

T, 


r-2 (T) = (AF|) =S p /T, (la) 

r a(T) = (^<*2) = S a /T . (lb) 

Here Sp and S ft represent the white noise strength for the variations in 
the X-ray flux F and the angular acceleration a. (The rms variations 
caused by fluctuations in the accretion process are denoted by 7p and T a , 
reserving a to represent other noise contributions such as those due to 
finite count rates.) Motivated by simple accretion models, we assume that 
concurrently sampled F-p and a<p are correlated, 

cov (F t , a T ) = prp (T) r a (T) , (2) 


with the correlation coefficient p close to unity, and that the correlation 
coefficient does not depend on the averaging time scale T. 
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Under ideal observing conditions of continuous sampling of F 
and a, the correlation between these variables can be estimated in a 
straightforward manner. The data are divided into a series of disjoint 
subintervals, each of length T, and concurrent averages and are 
computed for the kth subinterval. Mean values (F, a) and rms variations 
(rp , T a ) are estimated for the two variables F and a, and the sample 
correlation coefficient is computed using the statistic, 

f = ( F k-p) K-<*) , 3) 

The distribution of this statistic is relatively simple, giving a ready 
estimate of its uncertainty and providing a test for the significance of a 
nonzero value. 

A complication in applying this statistic to actual data is the 
presence of noise components in F and a in excess of the variations linked 
directly to the mass accretion process. We will denote the rms size of 
these excess noise components by op and a a , so that the total variances in 
the two variables are given by, 


var F = = 7p +Op , (4a) 

var a = = Ta+o# . (4b) 


To make the notation less cumbersome, explicit indication of the averaging 
time scale has been suppressed. 

As discussed by Boynton et al. (1986) and Boynton and Deeter 
(1986) , an excess noise component in the angular acceleration arises 
from fluctuations in pulse shape that induce white noise in pulse phase . 
For data sampled uniformly in time, the uncertainty in a from this noise 
component scales as T — , where T is the length of the analysis interval. 
On the other hand, the fluctuations in a arising from accretion torques 
are observed to have the character of white noise, for which the rms 
variation in a scales as T~ (Deeter et al. 1987, 1989). On short 
timescales, the uncertainty in a due to noise in pulse shape therefore 
dominates the torque variations. This situation is described by Boynton 
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et al. (1986) in terms of a crossover frequency in the power spectrum of 
variations in a, where the power law exponent is seen to change from zero 
to almost four — equivalent to the change of two in the exponent of the 
scaling laws for the rms variation described above. 

An important consequence of this crossover frequency (or crossover 
timescale) is that it restricts the averaging time appropriate for 
computing correlations. In particular, estimates of the angular 
acceleration on short time scales will be dominated by variations due 
to pulse-shape noise, and the correlation between F and a arising from 
mass accretion will be hopelessly obscured. This situation, wherein the 
excess noise increases rapidly at short time scales, can be compared to 
the case in which all sources of noise are approximately white and the 
rms variation scales as T -1 / 2 on all time scales. If this were true for 
both variables being correlated, then the correlation analysis could be 
undertaken on as short a timescale as practical, thus maximizing the number 
of correlation pairs. This freedom is lost in the present situation where 
the pulse-shape noise increases so rapidly with decreasing timescale, that 
it eventually dominates all other noise processes. 

Local values of the angular acceleration are obtained from analyzing 
the pulse phases, which constitute an available intermediate data record. 

We want to choose an method of estimating a which minimizes the ratio 
of the contribution from pulse shape noise to the contribution from 
fluctuations in the accretion torque. For example, in the case of 
constant acceleration, the contribution from pulse shape noise is minimized 
when the angular acceleration is estimated by the quadratic coefficient ctQ 
from a second-degree least-squares fit to the phase, 

<Kt) =<£o+fto(t-t 0 ) + ^(t-t o ) 2 , (5) 

where the coefficients <f> Q and 0 o (the phase and frequency at time t Q ) are 
also parameters in the least-squares solution. However, our model differs 
in that the actual acceleration is not constant but is characterized by 
white noise. In this case, the estimate of angular acceleration which 
achieves the greatest response to the actual acceleration with the least 


i 


( 12 ) 


response to pulse shape noise is the “lowest-mode estimator” discussed 
by Deeter and Boynton (1982). Using the methods of that paper, it can 
be shown that the estimate provided by the quadratic coefficient has an 
efficiency of 98% — that is, its variance is only 2% larger than that of 
the lowest-mode estimator. The determination of the lowest-mode estimator 
requires solving an eigenvalue problem, a procedure which is considerably 
more complex than the computation of the quadratic coefficient. In the 
interest of simplicity we therefore accept the latter as a practical 
alternative to the exact solution of the noise minimization problem. 

Having adopted the quadratic coefficient a Q as a simple and efficient 
local estimator of a, we next must assure that the estimates of F and 
a from an interval [a, b] are truly concurrent. This issue is treated 
in Appendix B, where we show that the estimate q q can be expressed as a 
sampling function g(t) applied to pulse phase. 


«o 




g(t)<£(t) dt , 


( 6 ) 


and that this integral can be rewritten using integration by parts to 
exhibit explicitly the equivalent sampling of angular acceleration, 


Q o - 
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The function g( — (t) is the particular second integral of g(t) which 
vanishes at the end points (a, b) , and is strictly continuous. On the other 
hand, g(t) can be highly discontinuous, jumping to zero wherever data is 
missing. Thus each local estimate of a from the quadratic coefficient is 
a continuous average of the entire run of a(t) in the interval, even if 
the flux and hence pulse phase are sampled discontinuously . 

The data set analyzed in this paper is quite discontinuous, as the 
the sampling is affected by Earth occultations, satellite operations, and 
other interruptions. Thus the average X-ray flux for a given interval is 
based on data from a subset consisting of many subintervals separated by 
significant gaps. The closest match for a flux average to the estimate 
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of a (eqn. 7) is obtained by applying the sampling function g — ^ (t) to 
the flux data on the subset of observations as indicated in Appendix C 
by equation (C3) . Thus the sampling of the two variables is not strictly 
concurrent; that is, there is a “mismatch” in sampling. In Appendix C we 
show that this mismatch leads to a dilution of the observed correlation 
coefficient with respect to the “true” correlation for concurrent sampling, 
even in the absence of excess noise in F and a. Indeed, any uncorrelated 
noise in either variable will further dilute the observed correlation. The 
diluted correlation coefficient p 1 is given by the expression, 


/>' = 


cov(F, a) 


= x'?^p, 


(var F var a)^-/^ r p r a 


( 8 ) 


where the factor x 1 accounts for dilution due to incomplete sampling of 
F (see eq. [Cl 1] in Appendix C) and the remaining factors account for 
dilution from excess noise (cf . eqs. [4a, b] , and eq. [C13] ). 


A final complication arises from irregularities in the data sampling 
that cause the dilution to differ among the various sample pairs (F^ak). 
To combine the correlation pairs within a single analysis requires that 
they be assigned weights in the computation. In Appendix D, we show that 
these weights should be proportional to the dilution yk = x k t 7 F/ T p] [ r ct / r ct] > 
and that a reasonable statistic for the correlation coefficient is 




k=l 


t' t> 

T f T a 


k=l 


The expectation and variance of r^ are 
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(10a) 

(10b) 


in the limit that the diluted correlation ykP is small for each summand 
(cf. eqs. [D8] and [D9] in Appendix D) . The variance of r^ is necessary 
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later on for determining the uncertainty in the estimated correlation 
coefficient and for evaluating the significance of a nonzero value. 

b~) Determination of Noise Strengths for F and a 

The determination of the white noise strengths Sp and S a (see 
eq. la,b) and the average values F and a have not yet been specified. 

The noise strengths correspond to the variances ?j£ and f& for the case 
of computing the correlation coefficient on homogeneous data (cf. eq. 2), 
and are usually computed directly from the actual scatter in the two 
variables . 

The case of X-ray flux is relatively straightforward, since 
fluctuations in addition to those contributed by the mass accretion 
process have similar noise properties. For the data we analyze, this 
additional noise component is negligible compared to the variations from 
mass accretion (cf. Appendix A), and an exact treatment is not necessary. 

We therefore absorb the additional noise component into the overall noise 
strength Sp . With this simplification, the variance for the kth interval 
scales inversely with the effective duration of the interval over which the 
flux is averaged, T^, 

r F,k - S F/T' k • (ID 

Thus a reasonable estimate for Sp is given by 

s F = ;E T k( p k- f ) 2 ’ d 2 ) 

k 

with a compatible estimate of F given by a weighted average of the 
individual F^ , 

f =XX F k/XX- (12) 

k k 

What complicates the computation in the case of the angular 
acceleration, a is the presence of a second noise source with distinctly 
different properties (white noise on <j> rather than a ) . The contributions 
to the variance of each estimate from pulse shape variations and from 
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fluctuations in mass accretion are comparable in magnitude for averages 
over intervals whose durations are close to crossover timescale, and so the 
pulse shape variations cannot be ignored. We could obtain an estimate of 
S a by subtracting the pulse-shape contribution to the variances of the a k . 
In practice, this “subtraction” is accomplished by choosing that value of 
S a for which the observed scatter just matches the expected scatter in 
the a k , 


£ 

k 


(a k - a)‘ 

v£k 



1 2 


k a,k 


(14) 


This particular method of weighting is obtained from the maximum likelihood 
criterion applied to S a .^ The average value of the angular acceleration a 


3 For a discussion of maximum likelihood estimates, see Lindgren 1962. 
The divisor T k and the extra factor in the denominators in equation 

(14) effectively reduce the importance of points with relatively large 
observational errors compared to intrinsic scatter. 


is specified by the equation 

k r a,k k r a.k 

Equations (14) and (15) constitute a non-linear problem in two unknowns, 
a and S a , that may be solved by a Newton- type method starting with 
S a = 0. Note that if the actual scatter is already less than the estimated 
observational errors, this procedure would give a negative value for S a , 
and therefore the contribution from intrinsic noise would be negligible. 

For each of the correlation sets presented here the observed scatter was in 
fact larger than that predicted from the observational uncertainties. 

In this section, we have shown that a mismatch in sampling between 
X-ray flux and pulse angular acceleration causes the measured correlation 
between them to be smaller in magnitude than the actual correlation. In 
addition, we have shown how to handle inhomogeneous data having different 
reductions for the various (F, a) pairs, by modifying the formula for 
the correlation coefficient to weight the pairs in proportion to their 
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dilutions. In any event, the reduction in the correlation coefficent must 
be taken into account in assessing the results of the correlation analysis. 


IV. CHOICE OF DATA SET FOR ANALYSIS 

The procedure outlined in § III for computing the correlation between 
X-ray flux and angular acceleration imposes some fairly severe restrictions 
on the character of data suitable for analysis. The X-ray flux and the 
angular acceleration should both have noticeable variations ascribable 
to the mass accretion process. Specifically, variations in the flux 
should be considerable larger than those due to counting statistics, and 
the X-ray flux should be sampled as completely as possible to minimize 
dilution of the observed correlation. Even though variations in angular 
acceleration due to mass accretion are expected to dominate variations 
induced by fluctuations in pulse shape on sufficiently long timescales, 
it is desirable that this timescale be as short as possible in order to 
obtain the largest possible number of correlation pairs. Finally, it is 
useful that the X-ray flux be observed in two or more energy channels, so 
that the effect of variable absorption can be assessed and removed. In 
particular, two energy channels, one above and one below a nominal 10 keV, 
are the minimal requirement to monitor variable absorption. 

Up to now, no data seem to have been obtained specifically for the 
purpose of studying F-ct correlations. Furthermore, most of the existing 
data sets on X-ray pulsars are clearly not suitable for this purpose. 

Early data sets, such as the extensive observations of Her X-l and Cen X-3 
by Uhuru, have low counting rates and uncertain collimator corrections, 
resulting in poorly determined X-ray fluxes and long crossover timescales . 
Later data sets from satellites with large detector area (such as HEAP ,1) 
were often limited to a day’s duration and were therefore too short to 
provide a useful number of correlation pairs. 

A promising X-ray source for studying F-a correlations from existing 
data is Vela X-l . Not only does this source exhibit large variations 
in both flux and angular acceleration, but it has also been observed 
extensively during the past ten years. These observations include 
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four fairly large sets: (1) 0S0 8 in 1978 May; (2) HEAP 1^ in 1978 

November — December; (3) Hakucho in 1980 — 1982; and (4) Tenma in 1983 
March and 1984 March. We will briefly discuss the merits and liabilities 
of these sets as concrete illustrations of the problems of designing or 
selecting a data set with favorable properties for an F-a study. For 
more detailed discussions of these data sets see Nagase et- al. (1984), 
Boynton et al. (1986) , and Deeter et al. (1987) . 

(1) The 1978 May set, obtained with QSQ 8, spans about 35 days and is 
characterized by a particularly high data recovery fraction of about 75%. 
However, a is well-determined only on intervals longer than 6 days, and 
after excluding eclipses there are only four correlation pairs available. 
Only one energy channel is available, and so there is no way of correcting 
for variable absorption. 

(2) The 1978 November — December data set ( HEAP consists of twelve 
half-day pointings over a span of 36 days, arranged in a hierarchical 
pattern with spacings from 1.5 days to 12 days. Determinations of 

a statistically different from zero are available on the relatively 
short interval of three days (coincidently the interval spanning three 
consecutive pointings at the 1.5-day spacing). Only two independent 
correlation pairs can be calculated on the three-day interval, with a 
data recovery rate of about 30%. Separate high and low energy fluxes are 
available . 

(3) The 1980 — 1982 set ( Hakucho) covers about 110 days in seven major 
blocks. Determinations of a are available on intervals longer than 
about 4 days, with a data recovery rate of ~40%. There are between 10 
and 20 independent correlation pairs, depending on the choice of interval 
length. Both high and low energy fluxes are available. This data set is 
our choice for analysis. 

(4) Observations with Tenma cover a total of 10 days in 1982 
March and 1983 March. Accurate determinations of a can be obtained for 
intervals longer than three days, with a data recovery rate of ~40%. 

There are three independent correlation pairs, using intervals of three 
days duration. Multichannel energy fluxes are available. 
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The single data set with the best properties is the Hakucho set, 
primarily because it affords the largest number of independent correlation 
pairs. As indicated above, we have to use intervals of length greater 
than four days in order to obtain sufficiently precise estimates of the 
a arising from fluctuations in rotation. Estimates of a on shorter 
intervals are dominated by pulse shape noise, thereby obscuring any actual 
correlation with the X-ray flux. On the other hand, the loss of nearly 

three days due to X-ray eclipse during each 8.9-day orbital cycle restricts 
the correlation analysis to intervals less than the six-day span between 
consecutive eclipses. Bridging eclipses causes unacceptable reductions in 
the data recovery rate, unless two consecutive orbits can be covered within 
a single interval. However, there are only four independent correlation 
pairs with the required interval length of 15 days, and this is too small a 
number on which to consider a correlation analysis. 

From these considerations, the time scales chosen for examination were 
6 and 4 days, represented by 11 and 17 correlation pairs, respectively. 


V. PRESENTATION OF COMPUTATIONS 

The observations of Vela X-l obtained with Hakucho have previously 
been studied by Nagase et al. (1982, 1984) and by Deeter et al. (1987). 
These investigations provide the initial steps in converting the raw 
counts to a form suitable for the present study of F-a correlations. In 
particular, we use the fluxes corrected for collimator shadowing determined 
by Nagase et aJL . (1984) , and the pulse phases computed by Deeter et al. 

(1987) . These reduced data require some further manipulations to covert 
them into F-a sample pairs. 

As explained at the end of § IV, we are constrained to test for 
correlations using an interval length between four and six days, and 
decided on the two extremes. The six-day intervals are entirely 
disjoint, but we allowed up to one day overlap for the four-day intervals, 
squeezing in wherever possible two intervals between consecutive X-ray 
eclipses. All of these intervals are listed in Table 1. 
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Although the fluxes and pulse phases were computed independently, they 
generally matched in time, with very few satellite orbits having data from 
one category and not the other. One criterion in selecting an interval 
for analysis was that both flux and phase data were both well distributed 
within the interval. In §§ Va and Vb we explain how these data were 
converted into estimates of the local flux and angular acceleration for 
each interval, and in § Vc we present the computation of the correlation 
coefficient and tests for significance of a nonzero value. 

a) Computation of Angular Acceleration 

The analysis of the Hakucho data by Deeter et al. (1987) provides 
residual pulse phases relative to chosen reference pulse frequencies. 

These phases have been corrected to the inertial frame of the mean orbit 
for Vela X-l reported in that paper. In principle, the residual phases 
could be converted to phases by adding back in the phase accumulated at 
the rate specified by the reference frequency. In practice, this step 
is unnecessary, since the addition or subtraction of a constant rate does 
not affect the second derivative (or quadratic term) associated with the 
angular acceleration. The average angular acceration a k for the kth 
interval was obtained by fitting the residual pulse phases A d> k j lying 
within the interval by a second-degree polynomial in time, 

Ad>(t) = Ad> k +Afi k (t - t k ) + ^a k (t - t k ) 2 , (16) 


with a separate fit performed on the phase data in each interval. The 
parameters A<^> k and Afi k are the residual pulse phase and frequency at 
the midtime t k (the value of a k determined from a second-degree fit is 
actually independent of the choice of t k ) . The parameters A(p k , AH k , 
and a are obtained by minimizing the sum of squared residuals between the 
observed and calculated phases in the k interval, 
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O 

weighting the phases by their inverse variances, o£ . . These variances 

k, j 

were determined from the scatter about local linear fits to the phase 
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data on intervals shorter than 1 . 6 days , and are reasonable estimates of 
the uncertainty induced from pulse shape noise (Deeter et al. 1987) . The 
least-squares solution which determines provides an estimate of the 
uncertainty in this quantity, also based on pulse shape noise. The values 
thus obtained for the and the associated errors are listed in Table 1. 

The solution to the least-squares problem provides an explicit 
expression for the coefficent in terms of the data j , 

a k = XX,j A< £k,j • ( 17 ) 

j 

For mathematical convenience, we transform this sum into an integral by 
introducing a function g(t) consisting of 6-functions, 

g(t) =XA,j 5 ( t - t k,j) • ( 18 ) 

j 

This function is used below in constructing the flux estimate F^ 
corresponding to . 

The mean a and the estimate of the noise strength S a (in excess of 
the noise attributable to pulse-shape noise) were then computed from the 
individual a^, using the method of § Illb, and the resulting values are 
reported in Table 2. 

b) Computation of Average X-ray Flux 

The Hakucho flux data comprises two energy channels, 1 — 9 keV and 
9 — 22 keV. There are obvious, strong variations in both channels at all 
time scales greater than a few hours (cf. Appendix A). These variations can 
be characterized as white noise with a high-frequency cutoff. In addition, 
both energy channels show systematic variations at the orbital frequency, 
which are stronger in the low energy channel. Detailed analysis indicates 
that these orbital variations are partly due to cold-matter absorption 
(Nagase et al . 1984) , and the remainder are presumably real luminosity 
variations. 
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For use in the F-e* analysis, the flux variations due to absorption 
should be removed. However, for the high energy channel, the absorption 
variations are less than 10% compared to total variations of more than a 
factor of five, and therefore need not be removed if only the high energy 
channel data are used in the analysis. There is, in addition, the 
possibility that systematic luminosity variations at the orbital frequency 
ought to be removed as well. This is because systematic variations in the 
angular acceleration at the orbital frequency (and its second harmonic) are 
removed in the orbital analysis, which relies on variations at just these 
frequencies to determine the orbital parameters (see e.g., Deeter, Boynton 
and Pravdo 1981) . The exact relationship between the orbital variations 
in the flux and angular acceleration are, however, model dependent. For 
instance, if the angular acceleration reverses sign sporadically (e.g., 
through accretion flow reversal), the positive and negative contributions 
to the variations in a at the orbital frequency would nearly cancel, 
while the corresponding variations in the flux would not cancel. 

With this warning, we first adopt the naive view of F-a correlations 
in which F and a vary linearly (or at least monotonically) , and the sign 
of a does not switch solely because of flow reversal. In this case, the 
orbital term in the angular acceleration due to mass accretion is removed 
by the analysis for orbital parameters. Thus we should also remove the 
systematic orbital variation in the flux. A simple analysis of the Hakucho 
data shows that this variation is closely in phase with the anomalistic 
cycle, with maximum near periapsis. To correct the fluxes in the high 
energy channel for this orbital variation, we multiply by the factor 

©(M) = 1 -0.25 cos M, (19) 

where M is the mean anomaly (i.e., measured from periapsis). 

We now turn to the question of computing average fluxes corresponding 
to the average angular accelerations determined in § IVa. In order to 
minimize dilution of the correlation between F and a, it is necessary that 
the sampling used to determine the average F on each interval match as 
closely as possible the sampling used in the determination of the average 
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a (Appendix C) . If we had continuous sampling of the X-ray flux, we would 
use the second integral g( — 2) (t) of g(t) for the flux sampling function. 

For the actual data, we have to take the gaps into account. 

The power density spectrum of the flux variations (Appendix A) 
suggests that there is a short time-scale cut-off in the flux variations at 
a few hours. We there adopt the procedure of counting any flux observation 
in a satellite orbit (duration about 1.5 hours) as providing full 
coverage for that orbit, assuming that the variation between consecutive 
observations is reasonably “smooth.” On the other hand, satellite orbits 
with no observations are treated as missing data. In the same vein, we 
consolidate all the data within a satellite orbit — usually covering 
about one-fourth of the 90 minutes satellite period — into a single flux 
average, correcting by the orbital factor 0(M) given in equation (19). 

The mean flux for the kth interval [a^b^] is then given by the 
discrete version of equation (C3) , 

**k = XI X^kJ )k(tk ? j )^k,j / XI x ^k,j )k(tk,j ) ’ (20) 

j j 

where the dummy index j ranges over all satellite orbits within the kth 
interval, t^ j is the mid-time of the jth orbit, F^ j is the average flux 
for the jth orbit, h(t) is the second integral of g(t), and x(b) is a 
“selector” function defined by 



for a flux observation in the jth orbit; 
for no flux observation in the jth orbit. 
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The uncertainty in the flux estimate due to photon statistics is very 
much smaller (at least a factor of 100 in variance) than the observed 
flux variations . In addition to uncertainties from photon statistics, 

it is likely that other variations in the observed flux which do not 
reflect actual variations in the luminosity. A primary cause for such flux 
variations (which was mentioned above) is variable absorbtion along the 
line of sight. These variations are again relatively small compared to the 
observed range of variations in the flux, but are difficult to quantify. 
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In Table 1 we give the uncertainties in the flux estimates Fj^ based solely 
on photon statistics. 

From the individual estimates of flux from each interval, we compute 
an over-all mean flux, and the strength of the white noise variations 
determined from the scatter in the sample {F^} (cf. § IHb) . The resulting 
values are reported in Table 2. 

c) Results of Correlation Analysis 

We computed correlation coefficients of the X-ray flux and angular 
acceleration at 6-day and 4-day time scales using equation (9) . These 
correlation coefficients (r^ ) are listed in Table 2, and are uncorrected 
for dilution. 

The uncertainty in the correlation coefficient is as important as 
the actual value. Also valuable is the distribution function for the 
statistic used to compute the correlation coefficient, since the test for 
significance of a non-zero coefficient requires knowing the distribution. 
Unlike the distribution for r (the statistic for a homogeneous collection 
of correlation pairs) which depends only on the actual coefficient p 
and the number of correlation pairs, the distribution for r* cannot be 
described in a simple fashion and so cannot be readily computed and 
tabulated. However, the statistic r 1 was designed to behave as closely as 
possible like the statistic r (with an “effective” number of degrees of 
freedom replacing the actual number of correlation pairs) , at least in the 
limit of near homogeneity for which all correlation pairs have comparable 
dilution. 

In the situation of homogeneous correlation pairs, the variable 
z = tanh p is approximately normal with variance approximately n — 3, 
where n is the number of correlation pairs (see, e.g., Bethea et al. 1984). 
We adopt this transformation for p^ , assuming near homogeneity. In this 
case, n should be replaced by the effective number of degrees of freedom v> 
and so the ±1 a errors for p ^ are given by 

Ar = tanh [tanh — ^ r 1 ± (y — 3) ^ ] — rt ? 
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with the usual convention of setting p t = r^ , as the most likely value 
for this parameter. These errors are listed in Table 2 along side the 
estimates for r ' . 

The undiluted correlation coefficient r / is readily computed from r^ 
by dividing by the mean dilution y (cf. eq. [D6] ) . Error limits on t 1 are 
obtained in a similar fashion. The resulting quantities are also reported 
in Table 2. 


As a check on the accuracy of the error limits derived with the 
tanh - * transformation, we also estimated directly the variances of r^ 
and t 1 by a scheme involving repeated resampling with replacement from 
each set of correlation pairs. Resampling schemes such as this are often 
called “bootstraps”; see, for instance, Efron (1982). An estimate of t 1 
is computed for each resampling, and the variance and standard error of 
the statistic t' for the entire sample is obtained from the scatter in the 
resampled estimates. The standard errors for T f obtained in this way are 
very close to the errors obtained from the distributional properties of the 
statistic t 1 (see Table 2) . 


In addition, we tested whether the hypothesis p = 0 can be rejected. 
For this test we again use the diluted statistic p t , whose statistical 
properties closely match those of the standard correlation coefficient. 
Under the hypothesis p* = p - 0, the statistic 


= r t 


-2 l 1 / 2 


rt2 
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has Student’s t-distribution with v — 2 degrees of freedom. The two-sided 
confidence level 0 is then determined by 


0 = Prob(|t| < |t obg | : p^ = 0) , 


(24) 


from a table of the t-distribution . 

The test for a correlation between flux and a is suggestive, but 
not compelling; the hypothesis that p = 0 is unlikely but cannot be rejected 
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at the 90% level for either the 4- or 6-day averages. However, the F-a 
correlation plots shown in Figure 3 have a suggestive “butterfly” pattern, 
with large values for F apparently correlated with both large positive 
and large negative excursions in a. We are thereby led to investigate 
F-|a| correlations, at the same 6-day and 4-day time scales. For these 
computations, we allowed the vertex of the butterfly to be located at a 
non-zero value of a, which is consistent with physical models. The 
value e*£ = —2.5 x 10 — ^ rad s~ ^ for the vertex was chosen to maximize 
the correlation between flux and |a — af|. We designate this case as “F-|a| 
correlations. ” 

In line with the discussion at the beginning of § Vb, this case most 
likely corresponds to the physical situation of flow reversal. In this 
case the observed orbital variation in the flux does not imply a net 
orbital variation in the angular acceleration which will bias the orbital 
parameters, since the orbital term alternates sign, depending on the 
direction of disc rotation. Thus the orbital flux correction was not 
applied, since in this case the corresponding variations in the angular 
acceleration presumably have not been removed in the orbital analysis. In 
the analysis for F-|a| correlations, the effective degrees of freedom were 
reduced by 2 to account for the decision to fold the angular acceleration 
and the determination of the folding parameter <*£ . With these minor 
chages, the F— jor| analysis follows that of the previous F-ct analysis. 
Numerical results are given in Table 2, and the resulting correlation 
diagrams are shown in Figure 3 along with the earlier results. 

The tests on the F-|a| correlations give a slightly more significant 
result for a non-zero correlation coefficent than do the tests based on 
F-a correlations. The most significant test is for the hypothesis p = 0 
applied to the F-|a| correlation for the 6-day intervals. This hypothesis 
is rejected at the 92% level, and this rises to the 96% level for the 
one-sided test, p < 0. This result has not been significantly weakened by 
performing a large number of distinct tests, and is very suggestive of a 
correlation between X-ray flux and rotational acceleration. 
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VI. DISCUSSION 


a) Flow Reversal 

The hint of a correlation between X-ray flux and the absolute value 
of the angular acceleration for Vela X-l suggests the possibility of flow 
reversal in the mass-accretion process for this system. We discuss here 
whether other observational evidence is consistent with flow reversal in 
this system, but a full treatment of the theoretical implications will be 
left for later papers (e.g., Lamb et al. 1988). 

For simplicity, we assume that the relationship between instaneous 
F and a is monotonic and close to linear for “normal” accretion flow (in 
the same direction as the neutron star rotation) , and that essentially 
the same relationship holds between F and —a for flow in the reverse 
direction. The form of the observed F-a relationship then depends largely 
on whether the reversal timescale is longer or shorter than the averaging 
timescale for a and F. If the reversal timescale is longer, then the F-a 
diagram will display two distinct lines corresponding to the two cases of 
direct and reverse accretion flow. However, if the reversal time scale is 
shorter, then the diagram will “fill in” (Lamb 1985) since each observed 
value for a will be an average over one or more reversals in the flow 
direction. The appearce of the F-a diagram for the analysis of Vela X-l 
presented here suggests that the reversal time scale is longer than the 
averaging time scale of 6 days. This is only a guess, and we do not have 
enough data (11 correlation pairs) to make a definite test for the size of 
the reveral timescale . 

Flow reversal will also have an affect on the power density spectrum 
of the fluctuations in the angular acceleration. We interpret the 
variations of the X-ray flux — almost 100% modulated on the timescale 
of several hours (see Appendix A) — as evidence for noise in the mass 
accretion process. Flow reversal causes an alternation in the direction of 
the angular acceleration of the pulsar, thereby introducing variations in 
the angular acceleration about as large as the variations from the noisy 
accretion process. The timescale for flow reversal appears to be an order 
of magnitude larger than the timescale for variations in mass accretion, 
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corresponding to a frequency bandwidth one-tenth as large. Although the 
power in the two noise components is comparable, the power density for the 
reversal fluctuations is ten times the power density for the variations in 
mass accretion, with an increase in power density at analysis frequencies 
shorter than f ss (27 r • 6 days) - * . The absence of a large jump in the power 
density spectrum of the variations in a is therefore evidence against flow 
reversal in Vela X-l. Unfortunately, we have no additional evidence at 
this time which will aid in reconciling these conflicting conclusions about 
flow reversal in Vela X-l. 

b) Guidelines for Future Studies 

Although the present investigation has not produced a decisive result 
regarding F-a correlations in Vela X-l, we can still use the experience 
gained to help establish guidelines for future F-a studies. 

Studies of the noise properties of both flux and angular acceleration 
should be made before an F-a correlation analysis is undertaken for a 
particular data set. Both variables should show significant variations 
whose statistical properties indicate that they could arise from the same 
basic process. In particular, the power spectra for F and a should be 
have similar behavior at sufficiently low frequencies, where variations in 
angular acceleration from mass accretion dominate variations induced from 
pulse shape noise . In the following general discussion, we assume that 
the power density spectrum of variations in a from mass accretion is flat 
(power-law with exponent equal to zero) , while the spectrum of variations 
from pulse-shape noise obeys a power law with exponent equal to four. 

An important concept in evaluating whether an F-a correlation 
analysis is appropriate for a particular data set is the “crossover 
timescale” T co , taken to be the time interval on which the rms variations 
in the angular acceleration from flucuations in mass accretion equal 
the rms variations induced by pulse-shape noise. We will define T co as 
the time interval on which an estimate of angular acceleration has equal 
variances from the two noise processes as given by equations (C4a) and 
(Cll) in Appendix C. In general, both variances scale with the interval 
length T: = S a K a /T^ (mass accretion fluctuations) and a# = S^K^/T 
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(pulse-shape noise) where K a and K ^ are constants whose values depend on 
the specific method of computing the angular acceleration. For an estimate 
of the angular acceleration obtained from the quadratic coefficient of a 
second degree polynomial fit to pulse phase, the timescale T co at which 
r a = a a i- s given by 

T co = (504) 1 / 4 (S < ^/S a ) 1 /4 ~ 4 . 74 ( s^/S a ) 1 / 4 . (25) 

This expression is closely related to equation (12) of Boynton et al. 

(1986) , which gives the “crossover frequency” at which the power densities 
of the two noise processes are equal, l/fco = 27r(S < ^/S a ) 1 / 4 . The small 
difference in the numerical coefficients in the expressions for T co and 
1 /f co reflects the distinction between the two estimates; f C o arises from 
a comparision of power densities at a precise frequency, while T C o comes 
from comparing power within a broad frequency band with l/T C o representing 
a typical frequency in that band. 

As indicated in § Ilia, the crossover timescale T C o sets a definite 
lower bound on the range of timescales appropriate for the study of F-a 
correlations. Moreover, the maximum feasible timescale is generally not 
much larger, since the number of independent correlation pairs diminishes 
as the analysis interval is increased. A compromise must be struck 
between either a large number of correlation pairs with a corresponding 
reduction in the expected correlation coefficient, or a small number of 
pairs with very little reduction. This tradeoff can be expressed in terms 
of the likelihood of detecting a correlation between F and a, when the 
correlation has a specified nonzero value p. This likelihood should be 
as large as possible, over the range of possible subdivisions of a given 
data set with total observing time T^ 0 t • This maximization problem is 
facilitated by applying the transformation z = tanh~^r to the correlation 
coefficient r; z is approximately normal with variance (n — 3)~^ (Bethea 
et al. 1985). We then maximize the ratio z/ (var z )V^, as the size of the 
analysis interval is varied. 

In implementing this maximization procedure, it is necessary to 
specify the true correlation coefficient p between F and a, as well as the 
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total data span . If the length of the analysis interval is denoted 

by T, then the number of independent correlation pairs is n = T^ 0 ^/T and 
the observed correlation coefficient is p/ [l+(T C o/T)^] » including 

the dilution due to pulse shape noise in a (cf. eq. [C12]). With these 
subsitutions in the ratio z/(varz)^/^, we obtain the expression 


l tot 


-3 


1/2 


tanh 


-1 


[1+(T co /T)4]1/2 ’ 
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which we maximize with respect to T. Unfortunately, the analytic 
maximization of equation (26) does not result in a simple expression 
in terms of the various parameters. However, in the limits of a 
small correlation coefficient (|p| less than 0.5) and a large number of 
correlation pairs (n larger than 25 or so) , an approximate solution can be 
obtained by replacing tanh - ^z by z and substituting n for n — 3 in varz. 
The resulting optimum timescale is T « 3^/4 t co « 1.3T co , to lowest order 
independent of both p and T^ 0 t . Calculation shows that the optimum length 
is shorter for small sample sizes (less than 25 or so) , and longer for 
large \p\. In any event the minimum is fairly broad, and any value of T 
between T co and 2T co will be close to optimum unless the sample size is 
very small or |p| is close to 1. 


This argument establishes the optimum length for the analysis 
interval, irrespective of whether a detection of an actual correlation 

is likely to be established. The minimum number of correlation pairs 

needed to establish a nonzero correlation clearly depends on the actual 
value of the correlation coefficent. To get a 2 a detection, we need to 
have z/(varz)*/^ = 2. Using the guideline that the analysis interval 
should be approximately 1.3T co and interpreting T^ 0 t/T as n, we can set the 
expression in equation (26) equal to 2 and solve for n, Thus, if p = 0.5, 

a minimum of 22 correlation pairs is required to assure a reasonable 

likelihood of a 2 o detection, with more or less pairs for smaller or larger 
p, respectively. 


The minimum data span required is approximately 1.3nT co . It 
is therefore essential to choose a source with a short crossover 
timescale, preferably a few hours, in order to keep the total data span 
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to a practical length. A small T co is dependent on a large (noise 
from the mass-accretion process) or a small (noise from pulse-shape 
variations) . 

Both conditions are more likely to be satisfied by a rapidly-spinning 
pulsar than by a slowly-spinning one. First, a rapidly-spinning pulsar is 
likely to be a “fast” pulsar in the sense of Ghosh and Lamb (1979b) , and 
hence is likely to have large torques and relatively large torque noise 
from the mass accretion process. Second, in considering contributions 
to S^, Boynton et al. (1986) have argued that the dominant pulse shape 
variations for high-count rate data is likely to be fluctuations intrinsic 
to the source as opposed to variations due to counting statistics. We have 
found that the size of pulse-to-pulse fluctuations in Her X-l and Vela X-l 
are nearly the same, even though the properties of these pulsar are quite 
different and their pulse frequencies differ by a factor of 250. If this 
property is found to hold for all pulsars, then the noise strength 
due to intrinsic pulse shape fluctuations would be inversely proportional 
to the number of pulses per unit time (i.e., the pulse frequency), and 
so should be smallest for the most rapidly-spinning pulsars. Thus the 
relatively rapid pulsars Her X-l , Cen X-3 and LMC X-l are more promising 
candidates for future study than the much slower Vela X-l. 


VII. CONCLUSIONS 

We have analyzed the usable Hakucho data on Vela X-l for correlations 
between X-ray flux and angular acceleration in the rotation of the neutron 
star. Our most significant result is a detection of a correlation between 
the flux and |a|, significant at the 9695 level. This result suggests the 
possibility of flow reversal in mass accretion for Vela X-l. We do not 
have any other definite, supporting evidence of this possibility, but the 
possibily of flow reversal in Vela X-l should encourage F-a studies of 
other X-ray pulsars which might also be fed by wind accretion. 

A major emphasis of this study has been the development of methods 
to handle peculiarities in the Vela X-l data which complicated a direct 
application of standard correlation analysis. The major difficulties 
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inherent in our data were a sampling mismatch between X-ray flux and 
pulsar angular acceleration; the dilution in correlation resulting from the 
mismatch; and the treatment of correlation pairs with different dilutions 
due to inhomogeneties in the data. The solutions to these difficulties are 
discussed with a fair amount of detail to allow use in future studies. 

The primary application of the methods developed here is to pulsars 
for which the angular acceleration can be determined only at time scales 
longer than the time scale for variations in the X-ray flux and angular 
acceleration. In some sources this may not be true — for instance, 
in the transient source EXO 2030+375 (Parmar et al. 1988) — and the 
statistical treatment presented here will not be necessary. (However, 
transient sources present other problems, such as ambiguities in 
separating accelerations in the pulse frequency due to orbital motion 
from accelerations in the rotation rate of the neutron star, particularly 
if the flare lasts for less than one binary orbit.) 

Finally, the methods discussed here are useful in planning future F-a 
studies, such as the amount of data needed to obtain an adequate number of 
correlation pairs and underlining the necessity for obtaining data with 
nearly continuous sampling of the X-ray flux to avoid a sampling mismatch 
with the angular acceleration. 
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APPENDIX A. 

POWER DENSITY SPECTRUM OF THE X-RAY FLUX FROM VELA X-l 


The fluxes used in the power-spectral analysis were from the Hakucho 
high-energy channel (9 — 22 keV) , normalized to center collimator. The 
high-energy channel was chosen over the low-energy channel (2 — 9 keV) 
because it was less affected by variable cold-matter absorbtion due to 
material in the vicinity of the source. The fluxes were given in counts 
per second, averaged over one pulsation period. These points were further 
consolidated to one point per satellite orbit, and data taken during 
eclipse were ignored. 

The power spectrum of these fluxes was computed using a variant of the 
simple periodogram: 


n n 

aj = ^ AFjj cos 27rftj c / ^ cos^ 27rftj c , (Ala) 

k=l k=l 


n n 

bf = ^ AFjj sin 2^ftj J ^ sin^ 27rftj c , (Alb) 

k=l k=l 

where the sums range over all data points in a particular interval. 

However, to avoid domination of the rest of the power spectrum by the 
relatively large constant term, we replace each F^ by its residual from the 
mean within the interval, 

AF k = F k - F , (A2) 


f = < A3 > 

k=l 

Likewise, we replace cos 2^‘ftj c and sin27rftj c by the residuals from each of 
their respective means, obtaining modified versions of equations (Ala) and 
(Alb) , 

n n 

a f = X> F k A cos 27r f t^ / ^ A cos^ 27rftj c , 
k=l k=l 


(A4a) 
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n n 

bf - E AF k A sin 27rftj c / A sin^ 27rftj c . (A4b) 

k=l k=l 


To convert the periodogram to a power density spectrum, we normalize 
by multiplying by the time covered by the observations, excluding essential 
gaps in the data record. In computing this time coverage, we also take 
into account that the flux variations are suppressed at time scales shorter 
than the minimum spacing of the data points At sa ^ = 0.067 days, the 
sidereal period of the Hakucho orbit) . Thus the total “live” time is given 
by AT = nAt sa t , and the power spectrum is 

P f = nAt(a|+b|) . (A5) 


The single periodogram determined from the entire span of data does 
not give an accurate representation of the power density spectrum at all 
analysis frequencies, since sampling periodicies transfer power between 
remote frequencies. In order to circumvent this difficulty, we divided 
the the data into subsets on a range of timescales with subset divisions 
chosen to provide reasonably uniform data spacing within each subset. The 
time scales chosen as nominal lengths for subsets were: 1000, 83, 20, 

10, 5, 0.5 and 0.25 days. Power estimates were computed at each of these 
periods, and for shorter periods in some cases F sae. TQ M-a==H 3r. There are 
two gaps in this hierarchy of time scales — around 360 days due to annual 
sampling of the source, and around 1 day due to daily sampling. 

All the power estimates at the same period were averaged, thus 
producing pooled estimated at 11 time scales. In addition, the noise 
contribution from photon statistics to these power estimates were computed 
from the corresponding errors, a ^ , on the individual flux estimates: 

n n 2 

var(a£ ) = <7^ cos 27rftj c / [ ^ cos^ 27rftj c , 

k=l k=l 


(A6a) 
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n n 2 

var (bj ) = ^ cr^ sin 2 7rf t^ / ^ sin^ 27rf , 

k=l k=l 

P f,phot = n At[var (a.f)+ var(bf)] . 

The resulting power spectrum is shown in the accompanying Figure 


(A6b) 

(A6c) 


In every case, the observed power density is much larger than that 
from photon noise. In addition, there is a deficiency in power density 
at the shortest periods, 1/f « 0.5, 0.25 days. This result confirms what 
is fairly evident in the light curves themselves — that the flux varies 
substantially (50% or more modulation) on the time scale of a couple of 
hours or longer, but the variations are much smoother on shorter time 
scales. [A formal fit with a Gaussian roll off gives a smoothing time 
scale of about 0.05 days.] There is no evidence in the power spectrum 
for variations around 5 — 10 days (around the orbital period) , but this 
is understandable because a narrow feature is greatly diluted in a 
broad-band version of the power spectrum. 


i 
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APPENDIX B. 

SAMPLING OF ANGULAR ACCELERATION 


An important consideration in testing for a correlation between flux 
and angular acceleration is selecting a method for estimating the angular 
acceleration from a given segment of data, and understanding how this 
estimate effectively averages the instaneous acceleration. A seemingly 
straightforward estimate is suggested by the low-order terms of a Taylor’s 
expansion, 

<f>(t) = 4 >o+Qq (t — to) (t — to)^ • (Bl) 

This expression indicates that the quadratic coefficent in a least-squares 
polynominal fit to the the phase is a reasonable estimate for the angular 
acceleration. Mechanically following this prescription is not enough 
for our purposes. We need to understand how the estimate of angular 
acceleration arises as an average of the instantaneous acceleration over 
the interval, in order to average the X-ray flux in a similar fashion. 
Furthermore, we want to insure that the estimate is efficient in the 
statistical sense; that is, that the variance contribution from pulse-shape 
noise (white noise in pulse phase) for this estimate of the angular 
acceleration is small compared to the pulse-shape variance incurred in 
alternative methods. 

A natural setting for discussing both of these points is an existing 
framework for estimating red noise in rotation (Deeter and Boynton 1982; 
Deeter 1984) , and for propagating red noise into parameter estimates 
(Appendix to Boynton et al. 1986). In the latter reference, several points 
of interest to the present discussion are examined. First, when data x(t) 
in an interval [a, b] are fit by a set of functions gj (t) with coefficients 

x(t) = (t)-H(t) . 

j 


(B2) 
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the solution which minimizes the sum squared residuals can be expressed as 
a vector of linear functionals applied to the data, 




gj (t)x(t) dt . 


(B3) 


Second, the functions gj (t) are “dual” to the fitted functions in the sense 
that the two sets of functions are cross-orthogonal, 


j gj (t)gk(t) dt = 6jk. 


(B4) 


Furthermore, if the first m functions are monomials in t up to degree m — 1, 
then as a special case of equation (B4) the remaining dual functions are 
orthogonal to the low-degree monomonials, 


J gj (t)t k dt = 0, 


for j>m+l, k < m — 1. 


(B5) 


Applying this result to the present situation, we consider fitting pulse 
phase data by a constant term, a linear term, and any third funtion g(t) , 


<t>(t) =4> 0 +fi 0 (t - t D )+^g(t). (B6) 

Then the parameter 0 is obtained from the data by using the dual function 
g(t) as a sampling function, 


0 = 


b 

g(t)4>(t) 


dt. 


Furthermore, g(t) is orthogonal to two lowest-degree monomials, 



g(t)t k 


dt = 0. 


for k = 0, 1 . 


(B7) 


(B8) 


But, as discussed by Deeter and Boynton (1982), this is exactly the 
condition needed to rewrite equation (B7) via integration by parts to show 
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explicity how angular acceleration [the second derivative of is 

sampled, 

/•b 

= (-1) 2 / g(- 2 ) (t)a(t) dt. (B9) 

J a 

The integrated terms drop out provided we choose the particular second 
integral which vanishes at the end points of the interval, 

g(~ 2) (a) =g("2) (b) =0. (BIO) 


We can now return to the original topic: how does the quadratic term 

in equation (Bl) sample the angular acceleration in the interval [a, b]? 

As an illustration, we take g(t) = (1/2) (t — t Q ) 2 where t Q = (a+b)/2, 
with continuous sampling of pulse phase in the interval. Then, by direct 
calculation we have 


g(t) 


15 ' 

4At^ - 


3 (t to ) 2 


- At 


2 


•> 


(Bll) 


and 


g(-2) 


(t) = 


15 

16At 5 L 


At 2 - (t-t 0 ) 2 ' 


(B12) 


where At = (b — a) /2 is the half-length of the interval. Note that 
g(~2) is a bell-shaped function within the interval [a, b] , which 
vanishes only at the endpoints. Using the coefficient of the quadratic 
term as an estimator of the angular acceleration is therefore an 
intuitively appealing weighted average of the instantaneous values of 
the acceleration over the interval. 


Furthermore, these basic virtues of g(“^) (t) carry over to cases 
of non-continuous sampling, as long as the sampling is approximately 
uniform. Even in the extreme case where the quadratic function is fit to 
a discrete set of points in the interval, the function g(“^) (t) will still 
be continuous, consisting of linear segments. The crucial implication 
of this result is that the instaneous angular acceleration at every point 
within the interval contributes to estimate of the average acceleration, 
even for discontinous sampling. 
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APPENDIX C. 

SAMPLING MISMATCH BETWEEN X-RAY FLUX AND ANGULAR ACCELERATION 


One major difficulty in the correlation analysis is obtaing averages 
for F and a which are strictly concurrent. A suitable average of the 
angular acceleration is the quadratic coefficient in a polynomial fit to 
pulse phase, and is effectively a weighted average of the instaneous a(t) 
over the entire interval even though discontinously sampled (Appendix B) . 
The instaneous F(t), on the other hand, may be used to compute an average 
flux that is valid only on the subinterval containing usable observations. 

We first consider the effect on the observed correlation from this 
mismatch in sampling, laying aside the consideration of effect of other 
sources of noise in either the flux or the angular acceleration. We 
consider the case where data is available only on a chopped-up subset 
of a specified interval [a, b] . This subset can be identified through a 
chacteristic function, 


X 



where the flux is observed; 
elsewhere . 


(Cl) 


The angular acceleration a is estimated by a quadratic function fit to the 
phases. As explained in Appendix B, this is tantamount to sampling a(t) 
continuously with a bell-shaped function h(t), 


a = 



a(t) dt . 


(C2) 


[h(t) is an appropriate second integral of a function g(t) which 
samples pulse phase discontinously.] For the purposes of computing 
and interpreting a correlation coefficient, it is necessary to match 
the sampling of F as closely as possible to that for a. Thus the same 
sampling function h(t) should be used for F(t), but by necessity it will be 
restricted to the subset of observations, 

F = J h(t)*(t)F(t) dt 


/ 
J a 


h(t)x(t) dt 


(C3) 



(39) 


The variances of these estimates for a and F can be computed from the 
sampling functions and the noise strengths, 




r 


Tp - Sp 


[h(t)] z dt, 


[h(t) X (t)] 2 dt 


[/. 


b l 2 

h(t)x(t) dt 


According to equations (la,b) in the text, the factors following 
can each be identified with an inverse time scale, 


(C4a) 

(C4b) 


S ex and Sp 





dt . 


(C5a) 

(C5b) 


In general, T ; < T < (b — a) . These inequalities are intuitively sensible — 
in particular, a is sampled with a bell-shaped function whose mean width 
is less than the length of the interval [a, b] , and the sampling of F is 
over an even shorter interval. 


The covariance between F and a is affected only by data within the 
subset of observations where they are both sampled. The easiest way to 
compute this covariance is to introduce averages a r and a 11 on and off the 
subset of observations, 


A 

J a 


h(t)x(t)a(t) dt 


f 

J a 


h (t) x (t) dt , 
■b 


h(t) [1 - x(t)]a(t) dt 


// 

' J a 


h(t) [1 — X (t) ] dt. 


(C6a) 

(C6b) 


Then a = xa^+Cl — x)a^, where x is the weighted fraction of time under 
observation, 

dt / [ h(t) dt . 

7 J a 


-/ 


h(t)x(t) 


(C7) 



(C8) 
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The variance for a 1 has the same scaling factor as that for F, vara ; = 
S a /T / . The covariance between F and a is then given by 


/i? \ /t? (SpS tt )^/^ 

cov(F, a) = cov(F, xa ) = xp — - — — f , 


where p is the correlation coefficient for concurrently sampled F and a 
(cf . eq. [2] in the text). The diluted correlation coefficient between the 
mismatched averages is 


Q 


* 


cov(F, a) 
T^ r a 


= x(T/T ') l l 2 p. 


(C9) 


As an approximation, T^/T ~ x (with equality when the data dropout is 
uniformly distributed on the interval), giving x(T/T , )^‘/^ ^ x^/^ as an 
approximation for the dilution factor. 


We now add the complication of variations in F and a beyond those due 
to fluctuations in the mass accretion rate, 


var F = rp +< 7 p , (ClOa) 

var ol = • (ClOb) 


Excess variations in the X-ray flux include noise due to finite count 
rate and unmodeled variations in detector area. These variations are 
likely to conform closely to white noise; for our data they are small and 
are ignored. On the other hand, the excess variations in a arising from 
fluctuations in pulse shape (including noise due to finite count rate) 
have a completely different character. These fluctuations are observed to 
be white in pulse phase (Deeter et al. 1987) , and so can be characterized 
by a noise strength which is distinct from . Indeed, the noise 
propagation law is different (cf. eq. [C4a] ) 


o 


2 

a 



dt . 


(CU) 


This, in fact, is the variance that arises in the least-squares 
determination of a in the quadratic fit to pulse phase. 
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The overall dilution, including the effect of excess noise, is given 
by an expression similar to equation (C9) with 7p and replaced by the 
total variances, 


P 


* 


cov(F, a) 

T 

1/2 

_2 1 
r F 

1/2 

T a 

(var F var a) 

T 7 . 


r 2 i ^2 

L r F^ a pJ 


_ T a+° a. 


1/2 


P - YP - 


(C12) 


In the case of homogeneous data, where both effects are the same for all 
samples included in the correlation analysis, the dilution factor is also 
the same for all samples. In this case, ordinary correlation analysis 
is applicable. However, our samples are sufficiently inhomogeneous to 
require a modification to incorporate weights into the formula for the 
correlation coefficent, as discussed in detail in Appendix D. 
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APPENDIX D. 

COMBINING INDEPENDENT CORRELATION PAIRS 


As discussed in the text, it may be possible to characterize 
correlated variations in X-ray flux and angular acceleration (arising 
from the mass-accretion process) by a correlation coefficient p and two 
noise strengths, Sp and S a . In this Appendix, we discuss modifications in 
a standard method for estimating p to accomodate data with less than ideal 
properties. 

For homogeneous data consisting of correlation pairs (F^a^) with 
identical statistical properties the natural statistic for investigating 
the correlation cofficient is 

r = (»/-) Sktrk-^H-a) m) 

rpr a ’ V ' 

where the index k ranges over the correlations pairs, and r a and 7p are 
the observed rms scatter in the two variables (see, for instance, Green 
and Margerison [1978]). The expectation of this statistic is the true 
correlation coefficent, E(r) =p (ignoring a bias of order 1/n which is less 
than the uncertainty in r) . As explained in the text, we have to modify 
equation (Dl) to accomodate several complications in our data. These 
complications include: 

(1) The correlated component is expected to have the character of 
white noise, and the variance of this component therefore scales 
inversely with the length of the averaging interval (eqs. [la,b] in 
the text). Due to sampling irregularities, the interval lengths in 
our sample are not identical, and so the individual variances are not 
uniform among the correlation pairs. 

(2) The observed correlation coefficient is diluted because of the 
mismatch in sampling between F and a (Appendix C) , and the dilution 
is different for the various pairs (F, a) . 

(3) There are excess variations in both F and a besides the 
fluctuations due to mass accretion, which cause a further dilution in 
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the observed correlation coefficient beyond that due to the mismatch 
in sampling. Particular care must be taken in accounting for the 
excess variations in the angular acceleration, since they have a 
distinctly different character from the intrinsic variations. 

Points (1) and (3) specify individual variances for each correlation 

pair, in contrast to the uniform variances Tp and indicated in the 

denominator of the right-hand side of equation (Dl) , and in addition point 

(3) specifies that the individual variances include a contribution from 

excess noise. We therefore introduce the total variances tL , and T f , , 

F,k a,k’ 

defined as 


J 2 _ _2 ,_2 

r F,k " r F,k +cr F,k ’ 

J2 _ 2 , 2 

T a.k ~ r a.k + Vk ’ 


(D2a) 

(D2b) 


where , and , represent the intrinsic noise due to mass accretion for 
r ,K Ct,K 

the kth subset of the data, and the corresponding er^’s the variance due 
to excess noise. The first modification to equation (Dl) is to replace 
division by rpr a outside the summation with division by ^ inside the 

summation. 


The first modification would suffice if all correlation pairs had 
a common correlation coefficent even if they had different variances. 
However, because of non-uniform dilution (points [2] and [3]) the observed 
correlation coefficent will differ among the pairs. To accomodate this 
complication, we further modify equation (Dl) to produce a weighted 
estimate of the correlation coefficient, 


r 


I 


k y k r a,k r F,k 


X>k. 

k 


(D3) 


where the divisors y^ correct for dilution in the various estimates, and 
the weights are to be determined by a minimum variance principal applied to 
r f . We therefore investigate the variances of the individual summands in 
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equation (D3) , 

, _ (a k - a) (F k - F) 

r k v J J 
y k r a,k r F,k 


(D4) 


For the purpose of deriving relative weights for the various r k we assume 
that the denominators yjc 7 ^ k r a k are known exactly; in this case, E(r k ) = p 
and 


var r' k = [1 - (y k p) 2 ]/y k ~ l/y| • (D5) 

The approximation in equation (D5) is valid provided the diluted 
correlation y k p is close to zero. The statistic t 1 with minimum variance 
is then specified by weighting the terms according to their inverse 
variances, w k = y k , thereby obtaining 



(a k - a) (F k - F) 
J J 

r a,k r F,k 



(D6) 


The expectation of this statistic is p and its variance is l/)P k y k , 
least in the limit of all quantities y k /> close to zero. 

There is a minor problem with the statistic t 1 which prevents it 
from being used directly as an approximation in statistical tests usually 
applied to r. This problem is illustrated by considering the limiting 
case of identical dilution y for all terms, whereby the statistic t 1 
reduces to r/y. The factor y is extraneous in the application of the 
usual statistical tests, which are based on the observed (i.e., diluted) 
statistic r. To alliviate this difficulty, we introduce the statistic r^ 
with a slightly different normalization, 


rt = Y, y k 

k 


(q k - a) (F k - F) 

r a,k r F,k 


X>k 

k 


The expectation of r 1 is given by 


(D7) 


E(r f ) = pJ2 v k/Y,Vk 

k k 


CDS) 
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wherein the factor y = can re g ar ^ed as the mean dilution 

for the statistic r'. Likewise, the variance is given by 

varrt = ^y2/[^ yk ] , (D9) 

k k 

in the limit that y p is close to zero. As with the statistic r, this 
variance is related inversely to the number of terms in the summation. For 
r i , however, only those terns with the least dilution carry full weight is 
the summation and the effective number of summands given by the inverse 
variance, !/= 1/var r^, is less than the actual number of correlation pairs, 
n. There is clearly better justification to use u rather than n for the 
number of correlation pairs when approximating the distribution for r^ by a 
distribution for r. 
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TABLE 1 

Flux and Angular Acceleration Data Used in Correlation Analysis 


^mid 

(JD- 

2,440,000) 


a Tp 

X 

/f 1 

i? a t ? b 

*x *x 

(days) 

(p rad s~ 2 ) (days) 

(flux 

coverage) 

(cnt s -1 cm -2 ) 


6-day Intervals 


4310.63 . . . 

4.32 

-1.9 ± 0.4 

1.28 

0.325 

0.100+0.001 

0.113+0.001 

4319.23 . . . 

4.47 

-1.5 

0.4 

1.12 

0.257 

0.075 

0.001 

0.074 

0.001 

4596.86 . . . 

4.30 

+0.7 

1.0 

0.56 

0.135 

0.107 

0.002 

0.112 

0.002 

4615.61 . . . 

4.90 

-4.4 

0.4 

1.59 

0.313 

0.122 

0.001 

0.138 

0.001 

4625.35 . . . 

3.15 

-4.9 

0.9 

0.95 

0.329 

0.084 

0.001 

0.101 

0.002 

4668.80 . . . 

3.84 

-2.5 

0.5 

0.92 

0.231 

0.134 

0.002 

0.144 

0.002 

4991.96 . . . 

4.41 

+3.4 

1.3 

0.87 

0.212 

0.186 

0.002 

0.221 

0.002 

5000.39 . . . 

3.78 

+2.8 

1.2 

0.90 

0.249 

0.166 

0.002 

0.163 

0.002 

5009.36 . . . 

3.75 

-0.4 

1.3 

1.11 

0.320 

0.066 

0.002 

0.069 

0.002 

5028.50 . . . 

3.80 

-2.3 

1.3 

2.14 

0.569 

0.078 

0.002 

0.089 

0.002 

5323.67 . . . 

3.74 

-2.4 ± 2.0 

0.25 

0.076 

0.066+0.004 

0.071+0.004 


4-day Intervals 


4309.98 . . 

3.03 

+1.6 ± 1.2 

0.81 

0.300 

0.108 +0.002 

0.119+0.002 

4311.61 . . 

3.13 

-5.9 

1.2 

1.15 

0.374 

0.094 

0.001 

0.112 

0.001 

4318.02 . . 

2.49 

+0.5 

1.7 

0.83 

0.358 

0.081 

0.002 

0.070 

0.001 

4320.62 . . 

3.27 

-0.5 

1.7 

0.59 

0.178 

0.082 

0.002 

0.097 

0.002 

4597.56 . . 

3.11 

-0.3 

1.5 

0.40 

0.140 

0.098 

0.003 

0.114 

0.003 

4614.35 . . 

3.04 

-2.4 

1.0 

0.74 

0.275 

0.113 

0.002 

0.104 

0.002 

4616.67 . . 

3.28 

-4.9 

0.7 

1.05 

0.362 

0.135 

0.002 

0.165 

0.002 

4625.35 . . 

3.15 

-4.9 

0.9 

0.95 

0.329 

0.084 

0.001 

0.101 

0.002 

4633.72 . . 

2.41 

+1.5 

2.6 

0.79 

0.370 

0.116 

0.002 

0.139 

0.002 

4668.10 . . 

3.07 

-2.0 

1.1 

0.76 

0.239 

0.129 

0.002 

0.125 

0.002 

4676.61 . . 

2.38 

-1.8 

1.3 

0.53 

0.218 

0.056 

0.002 

0.048 

0.002 

4955.13 . . 

3.12 

+0.8 

2.1 

0.90 

0.313 

0.145 

0.002 

0.150 

0.002 

4992.61 . . 

3.05 

+5.4 

2.3 

0.64 

0.228 

0.224 

0.002 

0.275 

0.002 

5000.43 . . 

3.73 

+2.5 

1.3 

0.88 

0.248 

0.165 

0.002 

0.163 

0.002 

5027.78 . . 

2.88 

-0.1 

2.8 

1.71 

0.618 

0.088 

0.003 

0.097 

0.003 

5029.24 . . 

2.64 

-4.3 

3.0 

1.34 

0.517 

0.063 

0.003 

0.077 

0.003 

5323.67 . . 

3.74 

-2.4 ± 2.0 

0.25 

0.076 

0.066+0.004 

0.071+0.004 


a X-ray flux without orbital correction. 

b X-ray flux corrected by multiplying by the factor 1 — 0125 cos M, where M is 
the orbital mean anomaly of Vela X-l. 
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TABLE 2 

Correlations between Flux and Acceleration 



F-e* Correlations 

F-|a| Correlations 


6-day 

4-day 

6-day 

4-day 

Quantity 

Intervals 

Intervals 

Intervals 

Intervals 

F (cnt cm - ^ s — ^) .... 

0.116 

0.119 

0.106 

0.118 

Sp (cnt cm~^ s— 1) . . . 

170 

162 

133 

114 

a (p rad s — ^) 

- 1.31 

- 1.25 

1.94 

2.45 

S a (10 — rad - ^ s — ^) . 

1.86 

1.76 

0.72 

0.26 

Number of (F, a) pairs . 

11 

17 

11 

17 

v (effective dof) .... 

10.4 

16.3 

8.2 e 

13 . 4 e 

rt 

0 40 "*”® 
U , 4 U — 0.34 

u,<3 — 0.26 

°' 58 - 0.11 

0 oq + 0.23 
u,; * y — 0.29 

y a 

r' (= r^/y) 

0.509 

°- 78 ±0.6 1 8 

0.485 

0 7 i ± 0.45 
- 0.54 

0.469 

1 26 "*"® 

1 •^— 0.77 

0.333 
i,A - 0.87 

a , h 

±0.86 

± 0.72 

± 0.63 

± 0.77 

t £ 

1.26 

1.37 

2.09 

1.55 

0.84 

0 d 

0.76 

0.80 

0.92 


Notes: 

a y is the average dilution of the corration due to sampling mis-match. 

k is the estimated error in r 1 from a ‘ ‘ jackknife * * computation (see 
text for details) . 

c t = r^[(t^ — 2) / ( 1 — r^)]^, approximate distribution given by Student’s 
t-distribution . 

^ Probability of rejecting the hypothesis p = 0. 

e The effective d.o.f. for the F-|a| correlations have been reduced by 2 to 
account for the decision to fold and for computing the optimum folding point. 
Without this reduction, i/ = 10.2 (6-day intervals) and 15.4 (4-day intervals). 


(48) 

REFERENCES 


Bethea, R. H. , Duran, B. S. , and Boullion, T. L. 1985, Statistical Methods 
for Engineers and Scientists. 2nd ed. (New York: M. Dekker) . 

Boynton, P. E. , and Deeter, J. E. 1986, in Proc. of the Workshop on Time 
Studies of X-ray Sources . ed. S. Hayakawa and F. Nagase (Dept, of 
Astrophyics, Nagoya Univ.), p. 1. 

Boynton, P. E. , Deeter, J. E. , Lamb, F. K., Zylstra, G. , Pravdo, S. H. , 
White, N. E. , Wood, K. S., and Yentis, D. J. 1984, Ap. J. (Letters) . 
283 , L53. 

Boynton, P.E. , Deeter, J. E. , Lamb, F. K. , and Zylstra, G. 1986, Ap. J. . 
307, 545. 

CRC Handbook of Tables for Probability and Statistics (Second Edition), ed. 
W. H. Beyer (Cleveland: The Chemical Rubber Company). 

Deeter, J. E. 1984, Ap. J. . 281^, 482. 

Deeter, J. E. , and Boynton P. E. 1982, Ap. J. . 261 , 337. 

Deeter, J. E. , Boynton, P. E. , Lamb, F. K. , and Zylstra, G. 1989, Ap. J. . 
336, 376. 

Deeter, J. E. , Boynton, P. E., Shibazaki, N. , Hayakawa, S., Nagase, F., and 
Sato, N. 1987, A. J., 93, 877. 

Efron, B. 1982, The Jackknife, the Bootstrap and Other Resampling Plans 
(Philadelphia: S.I.A.M.). 

Ghosh, P., and Lamb, F. K. 1978 Ap. J. (Letters) . 223, L83. 

1979a Ap . J., 232, 259. 

. 1979 b Ap_;_ jL, 234, 296. 



(49) 


Green, J. R. , and Margerison, D. 1978, Statistical Treatment of 

Experimental Data (Amsterdam: Elsevier Scientific Publishing 

Co. ) . 

Lamb, F. K. 1977, in Proc. 8th Texas Symposium on Relativistic Astrophysics 
( Ann. NY Acad. Sci. . 302, 482) . 

1985, in Proc. Japan-U.S. Seminar on Galactic and Extraealactic 
Compact X-ray Sources . ed. Y. Tanaka and W. H. G. Lewin (Tokyo: 

ISAS) . 

1989, in Timing Neutron Stars . ed. H. Ogelman and E. P. J. van 
den Heuvel (Dordrect: Kluwer Academic Press), in press. 

Lamb, F. K. , Pines, D. , and Shaham, J. 1978a, Ap. J. . 224, 969. 

1978b, Ap^ jL, 225, 582. 

Lamb, F. K. et al.1989, in preparation. 

Lindgren, B. W. 1962, Statistical Theory (New York: The Macmillian Co.). 

Nagase, F. et al. 1982, ISAS Research Note 178 (Tokyo: ISAS). 

Nagase, F. et al. 1984, Ap. J. . 280, 259. 

Parmar, A. N. , White, N. E. , Stella, L. , Izzo, C., and Ferri, P. (1988), 

Ap. J. . in press. 

Rappaport, S., and Joss, P. C. 1977, Nature . 266, 683. 

van der Klis, M. , and Bonnet-Bidaud, J. M. 1984, Astr. Ap. . 135, 155. 


(50) 

FIGURE CAPTIONS 


Fig. 1. Pulse phase and harmonic phases of the Vela X-l pulse, 
as a function of time, based on Hakucho data taken in 1981 January. The 
horizontal axis is Julian day, minus 2,440,000, and the vertical axis 
is pulse phase in seconds. The top panel shows the residual pulse phase 
with respect to a local, constant pulse frequency. The remaining panels, 
numbered (2) through (6) to indicate the individual harmonic numbers, show 
the residual harmonic phases with respect to a quadratic fit to the pulse 
phase. For reference, the pulse period of Vela X-l is 283 s. 

Fig. 2. (a) Power spectrum of the noise in X-ray flux for Vela-X-1, 

based on the high energy (9 — 22 keV) channel. The abscissa is the 
logarithm of analysis frequency in days - ^ , and the ordinate is logarithm 
of power density (cnts^ cm - ^ s — . (b) Expected power spectrum of the 

flux noise if it obeyed an f — ^ power law (white noise in the derivative of 
the flux) . 

Fig. 3. (a) Correlation between flux and angular acceleration for 
Vela X-l, averaged over intervals of 6 days duration. The abscissa is the 
X-ray flux (9-22 keV) , in counts cm ^ s 1 , corrected to center collimator, 
and also corrected for a term to remove orbital variation. The ordinate 

is the average angluar acceleration, in units of rad s . (b) Correlation 

between flux and |a — a c |, with flux and a averaged over the same intervals 
as in (a). The value a Q = —2.5 x 10 rad s * used in this figure was 
chosen to maximize the correlation between flux and |a — a Q |. The units are 
the same as for (a) , but here the flux has not been correction for orbital 
variation. (c) Same as (a), except that here the flux and angular (d) Same 
as (b) , here using the 4-day averages. 
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